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We study vortices in a radially inhomogeneous superfluid, as realized by a trapped degenerate Bose 
gas in a uniaxially symmetric potential. We show that, in contrast to a homogeneous superfluid, 
an off-axis vortex corresponds to an anisotropic superflow whose profile strongly depends on the 
distance to the trap axis. One consequence of this superflow anisotropy is vortex precession about 
the trap axis in the absence of an imposed rotation. In the complementary regime of a finite 
prescribed rotation, we compute the minimum-energy vortex density, showing that in the rapid- 
rotation limit it is extremely uniform, despite a strongly inhomogeneous (nearly) Thomas-Fermi 
condensate density ps{r). The weak radially-dependent contribution (oc V^lnps(r)) to the vortex 
distribution, that vanishes with the number of vortices Nv as -^ , arises from the interplay between 
vortex quantum discretness (namely their inability to faithfully support the imposed rigid-body 
rotation) and the inhomogeneous superfluid density. This leads to an enhancement of the vortex 
density at the center of a typical concave trap, a prediction that is in quantitative agreement 
with recent experiments. One striking consequence of the inhomogeneous vortex distribution is an 
azimuthally-directed, radially-shearing superflow. 



I. INTRODUCTION 

An essential defining property of a superfluid is the 
way it responds to an imposed rotation. Because a uni- 
form superfluid state can only support irrotational flow, 
one might expect a vanishing response to an applied ro- 
tation in a simply-connected domain. However, dating 
back to seminal works of Onsager and Feynmani? it has 
been understood that a superfluid can indeed rotate, but 
does so by nucleating localized topological defects (vor- 
tices) that carry discrete units of vorticity, thereby re- 
vealing the superfluid's macroscopic quantum nature. In 
a finite system and for a sufficiently low applied rota- 
tion rate H. < fici , the energetic cost of exciting even a 
single vortex is too high and the superfluid remains sta- 
tionary. For higher rates, vortices are nucleated into a 
regular hexagonal lattice that, on average, approximates 
a uniform response to the applied rigid-body rotation. 

In recent years, rapid progress in the field of confined 
degenerate Bose gases has led to the experimental real- 
ization of vortex lattices containing large numbers of vor- 
^j(,gg3,4j5,6,7,8j9 ^ Perhaps the most striking feature of these 
lattices is their apparent uniformity despite the strong 
spatial variation of the local superfluid density ps (r) im- 
posed by the trap. Despite some recent attempts^^, un- 
til very recentljiii no simple physical explanation for 
this uniformity (which has also been observed in simu- 
lationsiS) has appeared in the literature. 

To emphasize the apparent puzzle that such vortex uni- 
formity presents, one need not look further than the vor- 
tex state of a type-II superconductor. There, vortices are 
well-known to be strongly pinned by (i.e. attracted to) 
regions of suppressed superfluid density created by ma- 
terial imperfectionsi^. Based on this analogy, one might 
expect that vortices in a confined Bose gas would be re- 
pelled from the center of the trap (where superfluid den- 
sity and therefore vortex kinetic-energy cost is largest), 
and would congregate near the condensate edges, result- 



ing in a highly nonuniform and concave vortex distribu- 
tion. 

In this Paper we present a theory of vortices in a con- 
fined spatially-inhomogeneous rotating superfluid, pro- 
viding a simple physical explanation for and computing 
corrections to a uniform vortex array. We explicitly cal- 
culate the vortex spatial distribution fiy{r) in a trapped 
degenerate Bose gas characterized by Ps{r), showing that 
it is given by 



ny[r)~^— + c\/ Inps(r), 



(1) 



with c = (87r)~^ ln(e^^/^^Ci;), m the boson mass, f the 
vortex core size and w = mfl/h the scaled rotation veloc- 
ity. Since typically the spatial variation of the Thomas- 
Fermi (TF) condensate density (or Psir), see Ref. [Ij) 
takes place on the scale of the condensate radius R{^), 
the r-dependent correction in ny{r) (second term) to the 
uniform rigid-body result Uvo = mD./irh is subdominant 
in the thermodynamic and large rotation limits, vanish- 
ing as l/Ny with increasing number of vortices. More 
explicitly, for an isotropic harmonic trap, in the Thomas- 
Fermi approximation we find 



ny{r) ~ UyQ 



1 



K' 



27r (i?2 _ ^2)2 ^2^ 



1^7^, (2) 



where the condensate radius R{Vt) = Rq/ y/l — (fl/ilt)'^ 
is swelled by the centrifugal force beyond the Thomas- 
Fermi trap radius Rq, diverging as the rotational velocity 
approaches the trap frequency Qt ■ 

Hence we show that, indeed, consistent with 
experimenta^i^iSii^, the vortex density is highly uni- 
form and is well-approximated by the rigid-body re- 
sult. This uniformity can be most transparently un- 
derstood by ignoring vortex discreteness (valid at high 
rotation ratesi^^i^) and thinking in terms of the ener- 
getically optimum superfluid velocity Vs(r). In a nut- 



shell, for an arbitrary smoothly-varying superfluid den- 
sity Ps [r) , the superfluid velocity of the rigid-body form 
v^ = rjz X r, which corresponds to the uniform vortex 
density riuo — ^V x v^, always minimizes the boson's 
London free-energy. In terms of vortices, this nearly uni- 
form vortex distribution n^ ss n^o is a consequence of a 
balance between the spatial variation of the kinetic en- 
ergy per vortex and the vortex chemical potential, both 
of which scale with Ps{r). While it is energetically costlier 
to position vortices in a region where ps (r) is high (i.e. the 
center of the trap), the vortex chemical potential (con- 
trolled by ps{r)i^) is also high there, compensating and 
leading to an approximately uniform vortex density. The 
breakdown of the analogy with vortices in type-II super- 
conductors is therefore due to the difference in the spatial 
dependence of the vortex chemical potential in the two 
cases. While in superconductors this role is played by a 
uniform external magnetic field H, in trapped superflu- 
ids the vortex chemical potential is proportional to Ps (r) 
and thus spatially nonuniform. 

As we show here, a spatially-dependent correction to 
ny(r) in Eq. |^ arises from vortex discreteness and the 
related inability of the vortex state to locally reproduce 
the uniform vorticity v^ corresponding to rigid-body ro- 
tation. As shown rigorously long ago by Tkachenkc^^, in 
the case of a uniform superfluid the resulting energetic 
frustration cannot be reduced and the energetically opti- 
mum vortex distribution is a regular vortex lattice with 
density n^o- In contrast, in an inhomogeneous conden- 
sate the associated kinetic energy-density cost is radially- 
dependent and can be lowered by a nonuniform vortex 
distribution. As illustrated in Fig. ^ our analytical pre- 
diction for nvif), Eq- @, is indeed experimentally ob- 
servable and shows remarkable agreement with recent 
JILA experimental. Furthermore, as we show below, 
our prediction for the relation between nv{r) and Ps{r), 
Eq. (Q , can be more directly experimentally tested by in- 
troducing a known inhomogeneity to ps (r) (by modifying 
the trap potential) and studying the induced changes in 
ny(r). Despite the deviation of the vortex density from 
the rigid-body value, within the London regime, the feed- 
back effect on the condensate density Ps{r) is negligible, 
which remains well-approximated by the Thomas-Fermi 
expression^^'^^. 

One immediate consequence of our prediction of the 
deviation of fii,(r) from the rigid-body value n„o is a fi- 
nite superfluid flow in the reference frame of the lattice 
(rotating with frequency Q relative to the lab frame). 
Even more interestingly, the inhomogeneous vortex dis- 
tribution ny(r) implies that the resulting azimuthal su- 
perflow in fact exhibits a radial shear, as schematically 
illustrated in Fig. |21 

A prerequisite to the study of the thermodynamic vor- 
tex state under a finite imposed rotation (which was our 
primary initial goal) is the understanding of the single 
vortex problem. Although a number of interesting re- 
sults for the few-vortex problem have appeared in the 
literature -"■^'^'^^'^^i^'^i^^i^^, to our knowledge no explicit. 
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FIG. 1: (Color online) Plot of vortex density for the case of a 
harmonic trap as a function of radius, Eq. ||5J (normalized to 
riuo, dashed line), for the following fl and R values (labelled 
by the former): Q = .SGilt and -R = 49^m; f2 = .STilt and 
R = 31/^m; ri = .40nt and R = 25/im. Points denote data 
adapted from I. Coddington et al^. 




FIG. 2: Schematic plot of the mean superfluid velocity (ar- 
rows) along with vortex positions (circles) in the frame ro- 
tating counter-clockwise at il (in which the lattice is static). 
Thus, because Vs 7^ fiz x r (and has a magnitude less than it), 
the superfluid flows past the lattice in the clockwise direction. 



full solution for the velocity Vs(r) — {fi/m)W9{r) around 
a vortex in an inhomogeneous superfluid has been pub- 
lished. Studying this problem in the London approxima- 
tion (valid deep in the superfluid state), we find22i2S*2I, 
that the phase gradient at position r near a vortex lo- 
cated off-axis at position ro is given by 



V0{r) 
V0a{r) 



z X (r-rpj 

(r - ro)2 

z X Vp.,(ro) 

2/5s(ro) 



In 
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(3b) 



with V0a(r) giving the anisotropic correction to the stan- 
dard uniaxially-symmetric (homogeneous superfluid) re- 
sult (first term in Eq. \6a!fi ). This additional curl- free con- 
tribution to the superflow has interesting consequences 
for the dynamics of individual vortices in an inhomo- 
geneous superfluid. Recall first that, in the absence of 
external forces, a vortex moves with the local super- 



fluid velocity. For the case of a vortex in a spatially 
inhomogeneous superfluid, V6'a(r) makes an additional 
contribution to this superflow, leading to the remark- 
able result (also discussed for specific trap potentials in 
Rcfs. 21 , 22.23 24^ that a single vortex at radius r in an 
inhomogeneous superfluid will precess about the trap's 
axis of symmetry with frequency 



h \drP. 



2p. e 



(4) 



Such precession has been seen experimentally in Ref. |28| 
with a quality factor of order 10. We also demonstrate 
that a pair of vortices near the center of the trap exhibit 
center-of-charge rotation about the trap while also rotat- 
ing about each other. The precession rate reassuringly 
vanishes for uniform ps in the thermodynamic limit, in 
which even an off-axis vortex is a stationary state. 

In the opposite limit, for r far away from the off-axis 
vortex at Tq, the additional contribution to the superflow 
induced by condensate inhomogeneity has the form of a 
dipole field due to a opposite-charge vortex at ro and 
a vortex at the center of the trap. Hence this analyti- 
cal dipole contribution in the far field effectively shifts 
the vortex to the center of the trap, leading to a su- 
perflow that, far from the vortex, is axially symmetric, 
circulating around the trap center and described by the 
phase-gradient '^0 — ^4^. 

The rest of the paper is organized as follows. In SecHTI 
we use a mean-field model (appropriate for our interests 
here) for a rotating trapped Bose condensate to recall 
some standard results and to derive the London equations 
for the superfluid velocity in the vortex state. In Sec. lIIII 
we solve these equations to find the superflow around a 
single off-axis vortex in an inhomogeneous condensate for 
a number of experimentally-motivated condensate pro- 
files. We then use these results in Sec. II VI to derive pre- 
cessional dynamics in one- and two- vortex problems. In 
Sec. we study the many-vortex state and derive our 
primary result quoted in the Introduction, namely the 
relation between the vortex and superfluid densities. In 
Sec. IVII we derive a vortex lattice elastic energy in an 
inhomogeneous condensate and use it to calculate the 
Tkachenko mode equations (in the incompressible limit) 
for a spatially-varying condensate profile. We conclude 
in Sec. lVIll with a summary of our results and an outlook 
to future research. 



degree of uniaxial anisotropy which reduces the problem 
to two-dimensions perpendicular to the (z-) axis of ro- 
tation, although our work can be straightforwardly gen- 
eralized to three dimensions. The corresponding boson 
energy density, written in a reference frame rotating at 
frequency f2, is given by^S 



— |Vci>p + (l/(r)-M)|$P + f|$|^-fii. 



(5) 



where the rotation frequency fl plays the role of the 
"chemical potential" for angular momentum L^ with 



ni,^ = r2$^[-i?iz- (r X v)]$. 



(6) 



Here, ^(r) is the trapping potential, p is the boson 
number chemical potential, and g is the s-wave scatter- 
ing potential. Expressing the energy density in terms 
of the condensate density ps{r}^ and phase d{r), de- 
fined through the condensate field <I>(r) = y/pje'^, (which 
microscopically is the macroscopically-occupied single- 
particle wavefunction) we find 



■9„2 



[{v^r + p^ivey] + (i^(r) - p)p, + '-p: 



2m 
'ihfl{z X r) • [^/p^^^/p^ + ipsV6] 



(7) 



Deep within a dense Bose-condensed state, where /Os(r) 
is large, it is convenient to eliminate the superfluid den- 
sity from Eq. (Uj) by solving the corresponding Eulcr- 
Lagrange equation 



= --— V' 
2m'- 



-9p: 



3/2 



vpr-Vpr(ve)'] + (^(r)-A*)Vpi 

hfiy/^{i X r) • ve. (8) 



Neglecting derivatives in Ps{r) and replacing V0 by its 
approximate rigid-body value (see Sec. Ill Al below 1 V0 ~ 
ujz X r, we have the usual Thomas-Fermi (TF) resuHs^S 

Ps{r)^{p~V{r) + ^mn\')/g, (9) 

for the condensate density profile. For the experimentally 
most relevant case of a harmonic trap V{r) ~ ^mVijr^ 
with Q,t the trap frequency, Eq. ^ is uniaxially isotropic: 



II. MODEL OF A TRAPPED 
(INHOMOGENEOUS) ROTATING SUPERFLUID 

The starting point of our analysis is the ground-state 
energy density of an interacting rotating trapped Bose 
gas. To eliminate any explicit time dependence, it is 
convenient to work in the frame of reference in which the 
boundary conditions and the thermal cloud (the "nor- 
mal" fluid, playing the role of the proverbial bucket) are 
stationary'^. For simplicity, we focus on a trap with a high 



Ps{r)^{p~-m{n1-n'y)/g, (10) 

with the TF cloud radius R{Vl) (defined by the radius 
where Ps{r) vanishes) given by 



R{n) 



Ro 



v/i-(f^M)2' 



(11) 



and Ro = ^2p/inil^ the cloud radius at zero rotation. 
Using this solution of Eq. ||Sll inside Eq. {Tj) and dropping 



unimportant constant terms we find the London energy- 
functional for a rotating superfiuid 

T,^ f 

E=^ d\psir)[iVey ~ 2u;{i x r) • V^], (12) 
2m J 

with uj = mQ/h. 

Under an imposed rotation, a superfiuid turns by nu- 
cleating vortices. In the London limit these are point 
singularities at a set of positions r^ where the phase 9{r) 
is nonanalytic and satisfies 



V X V6{r) = 2Trny{r)z, 



(13a) 
(13b) 



with n^(r) the vortex density. Eqs. (|13a|l . (|13bp . together 
with the phase Euler-Lagrange equation (§§■ = 0, with 
E from Eq. lO) 



= V-{ps{r)V0)~ 
= V-ip4r)V9), 



LoVpsir) • (z X r) 



(14a) 
(14b) 



determine the superfiuid phase 9{r) and the correspond- 
ing superfiuid velocity Vs(r) — —'V9. Here, Eq. Ijl4b|) 
applies for the case of a static trap {oj — 0) or, for a ro- 
tating trap in the experimentally relevant case of uniax- 
ial trap symmetry (i.e. Ps(r) = ps{r)). It is important to 
note that the vortex positions Yi{t) (and therefore nt,(r)) 
are static in the frame rotating with the vortex lattice; 
thus, our time-independent expressions containing nt,(r) 
are defined in that frame. We shall assume that vortices 
are static in the frame rotating with the normal compo- 
nent, i.e., at frequency Q.. 
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3: (Color online) Magnitude of the phase gradient along 



a line intersecting vortices in a hexagonal rigid-body lattice 
(solid curve). For comparison, the dashed curve shows the 
rigid-body result. 




FIG. 4: Schematic depiction of the number of vortices A'^, as 
a function of the applied rotation rate SI of a finite superfiuid. 
The solid line is the expected actual number of vortices start- 
ing at the lower critical angular frequency f2ci whereas the 
dashed line is the rigid-body result (see main text). 



summarized by Eqs. H12l) . H13a|) and Ill4b|) . a nontrivial 
problem. For fast rotation (large fi, such that the co- 
herence length ^ is comparable to the vortex spacing a), 
the vortex state is dense and we can neglect vortex dis- 
creteness, approximating v<,(r) and rii,(r) by arbitrary 
smooth functions. Expressing E in Eq. (|12|l in terms of 
Vs(r) (and dropping a constant), we have 



£^-f / d\p,{r)(ws 



f7(zxr))2. 



(15) 



which is clearly minimized by the rigid-body solution 



riz X r, 



(16) 



corresponding to a uniform vortex density n„o = ^/tt = 
toO/ttTi. 

Clearly, this rigid-body result for Vs cannot be correct 
for an arbitrary rotation rate f2 for the simple reason 
that it corresponds to a uniform curl of V0(r), in con- 
trast to Eq. (|13a|l . This discrepancy can also be seen by 
comparing a numerically computed magnitude of the ex- 
act superfiuid velocity for a uniform lattice of vortices to 
that for the rigid-body result, Eq. H16|) . In Fig. Owe plot 
|V0(r)| satisfying Eq. p3a(l for a hexagonal array of vor- 
tices along a cut through several vortices, and compare 
it to the rigid-body phase gradient plotted as a dashed 
line. 



B. LoAver critical velocity of a rotating superfiuid 



A. Rigid-body result 

As noted in the Introduction, the vortex discreteness 
embodied in Eq. (|13b|) makes the minimization of E, 



Before embarking on a proper calculation of the vor- 
tex density, Uy (r) , that correctly incorporates vortex dis- 
creteness, we recall an elementary manifestation of the 
discrete nature of vortices in the opposite limit of low 



il, namely the existence of a lower-critical (imposed) ro- 
tation rate ri^i below which the superfluid "refuses" to 
rotateSl. To compute ilci , we simply compare the energy 
Eq. I|12|) of a single vortex at the origin to the system's en- 
ergy in the absence of vortices. The frequency where they 
cross is rici- For f^ > 0, the solution of Eqs. (|13a|l . (|14b|l 
for a vortex at the origin leads to a counter-clockwise 
circulation: 



V9= ^. 
r 



(17) 



For simplicity, we take for /Os(r) a simple form that is 
given by a constant po inside a radius R and zero outside 
this radius. Thus, inserting Eq. (|17|) into Eq. H12|l and 
evaluating the integral we have 



E 



2771 ^ ' 



(18) 



where ^ is the size of the vortex core derivable from 
Eq. (jHJ in the presence of a vortex and providing a short- 
distance cutoff for London theory. Since clearly for the 
case of no vortex E' = 0, we find a well-known result for 
a superfluid confined to radius R (see, e.g., Ref. 13^ 



flcl 



h R 



(19) 



which, unlike the analogous problem of type-II supercon- 
ductors (because of the absence of screening, i.e., an infi- 
nite London penetration length), vanishes in the thermo- 
dynamic limit. Since vortices do not appear for fi < flci, 
Eq. p9(l gives the first hint of the violation of the rigid- 
body result. This suggests that even for 17 > fici the 
actual number of vortices Ny{Q) is expected to be less 
than the rigid-body result to which it must asymptoti- 
cally approach in the large Q limit. Therefore, based on 
these general arguments we expect Ny{n) to follow the 
sohd curve in Fig. ^Si24 

Thus far, we have provided elementary calculations ex- 
hibiting the limiting behavior of a rotating confined Bose 
gas in the extreme largeiS and small fl regimes. Be- 
fore turning to the general many-vortex problem, in the 
next section, we first study the (prerequisite) single vor- 
tex problem in a spatially inhomogeneous superfluid. 



III. SINGLE VORTEX IN A SPATIALLY 
INHOMOGENEOUS SUPERFLUID 

In spite of its importance to understanding vortex 
states in inhomogeneous rotating condensates, the single 
vortex problem in an trapped rotating BEC has received 
relatively little attention, with most treatments using the 
standard homogeneous-condensate vortex solution as a 
starting point and focusing on dynamicgS2i^i'2SiSiS^4SSi2&. 
Two notable exceptions are works by Rubinstein and Pis- 
men^" and by Svidzinsky and Fettcr^^. However, because 
of the focus in Ref. |25| on the single- vortex dynamics. 



only asymptotics of the superfluid velocity near an off- 
axis vortex was worked out. On the other hand, although 
an exact solution for a single vortex was found in Ref. l2(l 
it corresponds to a very special (i.e. not motivated by any 
experimental realization) ps (r) . 

In contrast, our interest in computing the vortex den- 
sity distribution (Sec. El requires a full analysis of the 
superflow corresponding to a single off-axis vortex in an 
inhomogeneous condensate and this is the subject of this 
section. To determine the correct superfluid velocity, we 
study the Euler-Lagrange equation for 9 Eq. (|14bl) (which 
expresses the conservation of the boson number current) 
along with the vorticity constraint for a vortex at Tq: 



V X V9^2TT6'^^\r-ro). 



(20) 



For the most general problem, these equations must be 
supplemented by a boundary condition n- {ps(r)'V9)\fj — 
enforcing the vanishing of boson current through the 
boundary of the system. Consequently, the energy and 
dynamics of a vortex are strongly affected by both the 
condensate inhomogeneity and the nontrivial boundary 
conditions''^. However, for the case of a Bose gase 
trapped by a smooth potential (relevant to experiments), 
the condensate density Ps(r) vanishes at the cloud's 
boundary, automatically satisfying the above vanishing 
current condition. 

Before analyzing Eqs. Ijl4b|) and (|20|l in detail, we note 
that (inspired by the vortex solution V6' = ^,^_^ ^-^ for 
a uniform superfluid), an exact solution to Eq. Ijl4b|) is 
given by 



V9{r) 



Ps(ro)z X (r-rp) 
Ps{r) (r-ro)2 



(21) 



Although this solution fails to satisfy Eq. H20() (satisfying 
it only to leading order in spatial gradients of Ps(r)) and 
therefore is not a proper vortex solution, it does exhibit 
the qualitatively correct behavior that we shall verify in 
this section: V6' (and therefore the superfluid velocity) 
is generally larger where ps (r) is small and smaller where 
Ps(r) is large, relative to the solution for a uniform ps- 

For a full and systematic analysis, it is convenient to 
represent 



9ir) = 9,{v) + 9a{r), 
in terms of the known singular part 



9v{r) ~ tan 



y-yo 

X — Xq 



(22) 

(23a) 
(23b) 



corresponding to a vortex for a uniform superfluid den- 
sity, with 



V0„(r) 



z X (r - rp) 
(r-rp)2 ' 



(24) 



that ensures the vortex carries the correct topological 
charge in Eq. (|2nj). The analytic part 9a{r) (V x 



V0a(r) = 0) is chosen so that 9{r) satisfies the Euler- 
Lagrange equation Eq. Ijl4bp . which reduces to 



Vps ■ VOa + PsV^Oa = -Vps ■ VOy. 



(25) 



(Henceforth in this section we are restricting attention to 
the uniaxial symmetry case Ps(r) = Ps(^)-) A virtue of 
this approach is that it reduces the problem to the ana- 
lytical solution of Eq. (|25|l subject to a known "source" 
field — Vps • V6't,(r), Eq. (|24|l . In the case of a uniform 
Ps, Eq. lPf)|l is obviously solved by ^9a{v) = 0, with 
V0 = V9y{v) reducing to the well-known result Eq. (Pl|l . 
In the presence of a nonzero gradient of Psi^), V6'a(r) 
generally provides a nontrivial (but curl-free) correction 
to the superfluid velocity in Eq. H24|l . In addition, since 
for a smoothly- varying Ps{r) we expect that V0 is ap- 
proximately given by the uniform result Eq. (|24|) . this 
formulation is naturally set up for a systematic expan- 
sion in the small parameter V In p^ . 

Finally we note that for a vortex at the center of the 
trap (rg = 0), "Vps oc — f leads to a vanishing of the 
source term — Vps • V0i,(r) = 0. Consequently for an 
on-axis vortex V6'a(r) = and the solution reduces to a 
simple axially-symmetric result V0 = 'V9y = 22^ = Vc/). 

Thus, below we naturally focus on the nontrivial case 
of an off-trap-axis (ro 7^ 0) vortex. Although for a generic 
Psir) no closed form solution to Eq. (|25|l is available, 
a systematic asymptotic analysis in all relevant regimes 
determined by three length scales (the condensate size 
-R, the vortex position Tq and the displacement from the 
vortex r — Tq) is possible and is presented in Sees. IIII Al 
and IIII Bl (near a vortex |r — ro| ^ rp <C R) and in 
Sec. IIII CI f away from the vortex ro ^ |r — ro|). 
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FIG. 5: (Color online) Phase-gradient magnitude (solid curve) 
through a vortex away from the center of the trap, for the case 
of the Gaussian condensate profile Eq. I26II . The vortex is lo- 
cated at spatial position ro — (5, 0) with R — 5 characterizing 
the scale of the condensate. For comparison, the dashed curve 
depicts |V6'„| only, neglecting the analytic field. 



A. Gaussian condensate profile 

Before studying the general case of an arbitrary con- 
densate profile Psir), here we first analyze Eq. H25|l for 
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FIG. 6: (Color online) Phase-gradient magnitude (solid curve) 
through a vortex near the center of the trap, for the case of the 
Gaussian condensate profile Eq. I26II . The vortex is located 
at spatial position ro = (1,0) with R = 5 characterizing the 
scale of the condensate. For comparison, the dashed curve 
depicts |V0„| only, neglecting the analytic field. 



the simpler special case of a Gaussian condensate profile 



p,(r) = poexp(-rV2i?2), 



(26) 



with R characterizing the size of the condensate. For this 
choice, Eq. if^ reduces to 



-r • V0„ + R'V'Oa - r • Va„ 



(27) 



Focusing on the most interesting regime, we solve 
Eq. (|27|l in the vicinity of the vortex. Close to the vor- 
tex at ro, we can neglect the first term on the left side 
of Eq. I|27|l in comparison to the second, — r • V^a + 
R^W'^Oa — R^W^Oa, reducing the equation for 9a to 



V^ea = ^V9.-r. 
A solution to Eq. (|28|) is then easily found to be 



9a (r) 



2R^^' 



ro) • (ro X z) In ■ 



rol 



R 



(28) 



(29) 



where the denominator of the argument of the logarithm 
amounts to a choice of integration constant, fixed by 
matching to an "outer" solution. Our approximate choice 
above may be motivated physically by noting that we ex- 
pect 9a to exhibit spatial variation on a scale given by i?, 
as will be confirmed by more general calculations to fol- 
low. Although ro • (ro X z) = 0, we have written Eq. H29|) 
in this way to emphasize that close to the vortex, the r 
dependence in 9a arises through the r — ro combination 
only. 

To emphasize key qualitative features of the solution 
9{t), we reexpress the vortex solution in terms of the 
azimuthal angle between r — ro and ro. In terms of 
(j) the singular part of the solution is (recall Eq. ((230)) 
simply 9tj = (p. From Eq. I|29|l the analytic contribution 



0a (r) can be seen to be given by 0a(</') 
to 






csin (f), leading 



(30) 



with c= 2W'^o\y' 

An important (required) feature of this solution is 
that Oa{<j)) is indeed analytic, carrying vanishing winding 

/„ d(f>d,f,6a{(l)) — 0, and therefore preserving the quan- 
tization in units of 27r of the winding of the full vortex 
phase 9{(j)). The consequence of the nontrivial analytic 
contribution is that 9((j)) obeys this quantization condi- 
tion by winding more slowly near the center of the trap 
(i.e. where « tt) than away from it (i.e. where w 0) in 
agreement with our intuitive approximate vortex solution 
Eq. laj. 

As we illustrate in Fig. [3 the analytical contribution 9a 
has a small but qualitatively important effect on the su- 
perfluid velocity in the vicinity of a vortex. We compare 
the vortex superfluid velocity magnitude | V0| for a vor- 
tex in an inhomogeneous condensate to its counterpart 
\'V9y\ in the homogeneous condensate along x through 
the vortex located at rg = xqX. While both velocity pro- 
files exhibit the characteristic l/r divergence, consistent 
with our general arguments the superfluid velocity |V0| 
is clearly smaller than \V9y\ near the center of the trap 
and larger than it away from the center of the trap, with 
this effect vanishing (see Fig. O for a vortex near the 
center of the trap. 



B. Generic spatially-varying condensate profile: 
"inner" solution 

We now turn to the full problem of finding the su- 
perfluid velocity near a vortex in an inhomogeneous su- 
perfluid deflned by a general (uniaxially symmetric and 
smooth, with i? 3> ^) Psir)- 

Our method is a simple generalization of the calcula- 
tion for the Gaussian case in Sec. IIII Al and applies in 
the regime |r — ro| <C tq. Working to lowest order in the 
gradient of the condensate profile, near a particular vor- 
tex we approximate the smoothly varying Ps{r) and its 
gradient 'Vps(r) by their values at the vortex position Tq. 
This considerably simplifies the Euler-Lagrange equation 
(while retaining the essential physics) to 

VM(ro) • V0„(r) + iv20,(r) = -Vp{ro) ■ V0„, (31) 



where we defined 



Vfi{r) 



1 Vps(0 

2 Ps{r) ■ 



(32) 



The solution to Eq. H31|) can now be easily obtained by 
first finding a corresponding Green function that satisfies 



(VMro)-V + -V2)G(r 



r')=(5(2)(r-r'). (33) 



An explicit expression for G'(r) (verifiable via direct sub- 
stitution) may be expressed in terms of the Bessel func- 
tion Ko{x): 

G(r) = -ie— ^^('^«)iro(r|V/i(ro)|). (34) 

TT 

This then leads to the solution to Eq. (|31|l given by 



a(r) -- / d2/G(r - r')V/i(ro) 



z X (r' - rp) 
(r' - ro)2 



, (35a) 



where in going to Eq. Ij35b|l we have shifted the integra- 
tion variable r' ^ r' -I- r, with the resulting integral on 
the right side clearly a function of r — Tq only (a con- 
sequence of our approximation above). Moreover, since 
the second factor is sharply peaked near r' ~ rg — r, it 
is valid to this order of approximation to replace G(— r') 
by its value at this point. The remaining integral may 
be easily evaluated, finally yielding 

9a{r) « -(r-ro)-(zx VAi(ro))e-('-"-°)-^^('-«) 

xKoi\r-ro\\V^i{ro\), (36a) 

« (r-ro)-(zx VM(ro))ln|r-ro||VM(ro)|,(36b) 



where in Eq. Ij36b|) we have taken the r ^ Tq limit of 
Eq. (|36a|l . Since p,{r) varies on the scale of the condensate 
size R, the analytical correction 9a scales like l/R^, and, 
as expected, vanishes in the thermodynamic limit. 

Using V^(ro) — —ro\drp,{ro)\ we find that at fixed 
distance near the vortex, the azimuthal (j) dependence of 
9{(J3) is given by 



0(0) 



COS (p 



(37) 



with c = \drP-iro)\\r - roliiTodr - i"o||9r/i(^o)|) and c' = 
|r — ro||9rM(''o)| positive ^-independent functions that in- 
crease with increasing |r — ro| (in the stated small |r — ro| 
regime). Because of the smallness of c' near rg, for a 
Gaussian condensate profile this general result reduces 
to that found in Sec IIII Al Again, as required by analyt- 

icity of 9 air), it is easy to show that /p d(l)d^9a[(j)) — 0, 
thereby preserving quantization of circulation of super- 
flow encoded in /^ d(j)d^9{(f>) — 2tt. 

Although the phase (and corresponding superflow) dis- 
tortion given by 9{(j)) is quite small (especially near the 
vortex) we note that there are a number of experimental 
techniques that have been used succesfuUy to measure 
the condensate phase growth around a vortex'^'"^^, giving 
hope that the prediction in Eq. H37|l may be experimen- 
tally testable. 

Using Eq. (|36b|) . the superfluid velocity near a vortex 
(defined by |r — ro||V^(ro)| ^ 1) is easily computed, 
giving 



v.(r) 



?i rz X (r 
m 



ro) 



L (r-ro)2 
zxVM(ro)ln(|r-ro||VM'^o)l) ,(38) 
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in agreement with the result found in Ref. ^. Hence 
quite generaUy, near a vortex the correction to the su- 
perfluid velocity arising from the inhomogeneity of the 
condensate is approximately spatially uniform and is or- 
thogonal to the displacement ro from the trap center. 
We see that, as anticipated based on general arguments 
above, the superflow is slower at (/> = tt and faster at 
(/) = (see Fig. [7|. Moreover, for Tq — > 0, V^(ro) -^ 0, 
so that the analytic field is negligible for a vortex near 
the center of the trap. 



/ 



\ 



/ 



FIG. 7; Schematic picture of the effect of the analytic field 
on the superfiuid velocity of an off-axis vortex (Eq. 1381 '): the 
vortex location is denoted by a filled circle. The solid arrows 
denote V6 as affected by the nearly uniform analytic field. 
On the innermost ring the nearly uniform V^a is depicted as 
a small dotted arrow and V^i, is depicted as a larger dashed 
arrow; their vector sum leads to the distorted VS that is 
smaller in magnitude near the center of the trap and larger 
near the edge. 



Generic spatially-varying condensate profile: 
"outer" solution 



We now calculate the asymptotics of the superfiuid ve- 
locity far away from the vortex. Interestingly, after some 
examination of the form of Eq. (|25|l it is clear that the 
outer solution can be easily found. To see this first note 
that for 6a that satisfies V'^Oa = Eq. (|^ reduces to 



Vps ■ V0a = -Vps ■ Vt 



(39) 



Naively this would be satisfied by the antivortex (at rg) 
solution 



__ z X (r-rp) 
(r - ro)2 



(40a) 
(40b) 



(which satisfies the assumed condition V^^a = 0) were 
it not for the analyticity requirement <f dv ■ "VOair) — 0. 
However, this is easily fixed by adding to the solution 
in Eq. Ij40bp another contribution ^4^ due to a positive 
vortex located at the origin, which clearly satisfies the 



homogeneous part of Eq. (|25|l . We thus obtain an ex- 
act outer solution"^^ (verifiable by direct substitution) to 
Eq. 123): 



V0a(r) 



ro 



(r-ro) 



(41) 



Because the solution Eq. 1411) corresponds to a vortex 
dipole, with a positive vortex at the trap center and a 
negative vortex at the location of the true vortex Tq it 
is easy to see that, as required, it satisfies the vanishing 
circulation (V x V0q = 0) condition for |r — ro| ^ tq. 
Hence far from an off-axis vortex at Tq, the superfiuid 
velocity is given by 



Vs(r) 



h z X r 



(42) 



i.e., the superfiow adjusts to be axially symmetric about 
the trap center, such that a vortex at r = Tq appears to 
be sitting at r = 0. 



2-K V 




FIG. 8: Plot of the phase around a vortex in an inhomoge- 
neous superfiuid for the Bessel profile Eq. 1431 as a function 
of the angle 4> between r — ro and ro, for Eq. 1441 integrated 
numerically. The parameters _R = 10 and ro = 20 and the 
curves represent |r — ro| = 1 (solid), 10 (dashed) and 100 
(dot-dashed). 



D. Bessel condensate profile: Exact solution due to 
Rubinstein and Pismen 



Although we have been unable to find a closed form 
vortex solution for an arbitrary condensate profile ps (r) , 
for a specially contrived but qualitatively realistic profile 



Ps{r) 



Po 



h{r/Rf 



(43) 



a closed form solution was found by Rubinstein and Pis- 
men in the context of optical vortices in an inhomoge- 
neous laser beam.^°. In Eq. (|43|l . Iq{x) is the modified 



9(0) - 0(0) 



2-Rl 
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FIG. 9: Same as Fig.|Hlbut with R = 20, ro = 5 and |r-roj = 
1 (solid line) and |r — ro| = 1000 (dashed line) exhibiting the 
smallness of the analytic field for vortices close to the center 
of the trap. 
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FIG. 10: (Color online) Contour plot of constant phase gra- 
dient curves for a vortex at ro = (1,0) with R = 10. Their 
near-circular symmetry indicates the unimportance of the an- 
alytic field in this regime. 
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FIG. 11: (Color online) Contour plot of constant phase gradi- 
ent curves for a vortex at ro = (10, 0) with R = 10. The effect 
of the analytic field is clearly seen in the distorted shapes of 
these curves; the phase gradient has been decreased to the 
left of the vortex and increased to the right. 
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FIG. 12: (Color online) Phase-gradient magnitude (Eq. I44II ') 
through a vortex far from the center of the trap, located at 
spatial position ro = (10,0) with R — 10. For comparison, 
the dashed curve depicts |V0^| only, neglecting the analytic 
field (solid curve). 



Bessel function, with the asymptotic behavior /o(a;) — *■ 
1 -|- x'^/4: and exp(x)/v27rx for a; ^ and x -^ oo, re- 
spectively. This asymptotics shows that this special ps (r) 
behaves as 1 — r^ /2B? for r ^ i? (i.e. it is parabolic) and 
decays exponentially for r ^ R. Thus, despite its fine- 
tuned very specific form, it is quahtatively similar to the 
experimentaUy relevant TF profile. 

Although most of the properties of a vortex are con- 
tained in the asymptotic solutions found in previous sec- 
tions, it is useful to check the predictions found there for 



a general ps (r) against the specific but exact solution for 
Ps{r) in Eq. H43|l . This is the subject of present subsec- 
tion. 

Following Ref. |23, the Euler-Lagrange equation (J14b|) 
can be solved exactly by introducing an auxilary function 
X related to 9 by ps{r)^9 = z x Vx. With this transfor- 
mation and the very special choice of ps (r) in Eq. H43|) 
the vortex quantization condition Eq. (|20|l reduces to 
a simple Helmholtz-like partial differential equation for 



10 



-1/2 

Ps x(^)- which in turn leads to the solutioi 



.20 



ve{r) 



^ii{ro)-l-i(r) 



ix [A'o(|r-ro|/i?)V/i(r) 



+Vi^o(|r-ro|/i?)], 



(44) 



where iJL{r) = — ln/o(r/i?) and Kf){x) is the Bessel func- 
tion. In the regime r,rQ <^ R of primary interest to us 
we find (using Kq{x) ~ — Inx, for a; -^ 0) 



ve{r) 



In- 



R 



z X f Ii{r/R) 
R loiro/R) |r-ro| 
i: X (r-rp) Io{r/R) 
(r-ro)2 loiro/R)- 



(45) 



For r ^ rg, the second term dominates inside V0, leav- 
ing the first term as a small correction. In this regime, 
the first curl-free term is the V0a analytic contribution 
introduced in Eq. (|22|l . Thus, this exact result agrees 
with the general one found in Eq. (|38|1 and discussed in 
the Introduction, Eq. H3a|l . 

In the opposite regime, far away from the vortex (r ^ 
rg), but still away from the condensate boundary (r <C 
R), the first term continues to be subdominant to the 
"bare" vortex contribution (second term). However, its 
r dependence can no longer be ignored and leads to a 
superfiow that is perpendicular to the position vector r, 
rather than the vortex position vector Tq. 

Even further in the far field, defined by r 3> i? but 
ro <C i? (so that the vortex is still near the center of the 
trap)22i, the exact solution Eq. ll^ reduces simply to 



V6'(r) 



z X r 



(46) 



with small power-law corrections in R/r, in agreement 
with our "outer" solution Eq. (|41|l . 

The vortex distortion due to the condensate inhomo- 
geneity can also be seen in the angular coordinate depe- 
dence of 6{(j)) — 9{0) (obtained by numerically integrating 
Eq. H44I) along a contour between and (p) that we plot 
in Figs. IHl and ini for different radii |r — ro|. As we found 
in Sec. IIII Al the distinction between the exact vortex 
6{(t)) and Oy{ip) — cj) is largest for vortices far from the 
center of the trap (see Fig. ISJ, and grows with |r — ro|. 
For a vortex near the center of the trap (Fig.^, the dis- 
tinction is barely noticeable. From examining these plots 
it is clear that an off-axis vortex exhibits an azimuthal 
dependence of the phase that (as found in Sec. IIII B|l is 
given by 9{(j)) « -I- c sin with c a positive constant that 
grows with distance away from the vortex, |r — ro|, and 
with distance of the vortex from the trap axis, vq. 

To further depict the increasing importance of the an- 
alytic field for vortices located away from the center of 
the trap, in Figs. 1101 and 1111 we plot (for R — 10) con- 
tours of constant |V6'| for a vortex at Tq — (1,0) and 
ro = (10,0), respectively. As expected from earlier dis- 
cussion, Eq. (|38|l , the superfiow of a vortex near the cen- 
ter of the trap (Fig. I10|l characterized by nearly circu- 
lar contours contrasts strongly with that of a far-off-axis 



vortex (Fig. 111(1 described by highly distorted contours 
of contant superfiuid velocity. This latter distortion is a 
manifestation of our earlier finding of a larger correction 
to the superfiuid velocity for vortices away from the cen- 
ter of the trap and smaller near it. In Fig. ^]wc depict 
the phase-gradient magnitude (solid curve) for the same 
parameters as Fig. ^2 contrasting it with that of a vor- 
tex in a homogeneous condensate (dashed curve) given 
by \VU 



IV. PRECESSIONAL DYNAMICS OF 

INDIVIDUAL VORTICES IN SPATIALLY 

INHOMOGENEOUS SUPERFLUIDS 



Having found in Sec. IIIII the vortex superfiow profile 
in a spatially inhomogeneous superfiuid, here we study 
implications of these results for few-vortex dynamics. In 
the absence of rotation (17 = < Jlci), vortex excita- 
tions are metastable, and therefore should eventually es- 
cape the condensate, allowing the superfiuid to relax to 
the ground state. However, because of the (massless) 
guiding-center dynamics of the vortex, which conserves 
angular momentum, we expect the vortex lifetime to be 
quite long (limited by the trap's azimuthal asymmetry 
and decoherence) . As we show in the next subsection, 
on shorter timescales a vortex in a confined Bose gas, 
located away from the trap axis, necessarily precesses 
around the center of the trap22iS2i^. This can be under- 
stood from general considerations by noting that energy 
eigenstates in an axially-symmetric (finite) trap are also 
angular momentum eigenstates. Consequently, an off- 
axis point vortex (that clearly is not an eigenstate of an- 
gular momentum) is not an energy eigenstate and must 
therefore evolve at constant energy. As we demonstrate 
in Appendix ^ such vortex precession is a property of 
any confined superfiuid, even with a uniform condensate 
density. 



A. Precession of a single vortex 

We shall now demonstrate (via two different routes) 
that an off-axis vortex in a spatially inhomogeneous su- 
perfiuid precesses around the trap centeri22i22i24. To be- 
gin, we recall that the dynamics of a vortex is governed 
by the Magnus "force"i224a 



'M 



2TThpsi X {Vi ~ Vs 



(47) 



where r^ is the velocity of the ith vortex and v^ is the 
background superfiuid velocity (not including the sym- 
metric divergent part V^i, due to the vortex itself) at 
the location of the vortex, r^. In the absence of other 
forces and in equilibrium Fm — 0, giving that a vortex 
moves at the local superfiuid velocity, r^ = v^. In an 
inhomogeneous superfiuid, it is appropriate to identify 
Vs with the curl-free contribution K^Oa/Tn that, as dis- 
cusssed in Sec. IIIII arises from condensate inhomogeneity. 
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Thus, using the second (analytic) term of Eq. (|38|l with 
|r — r^l « ^, we find 



—zxV^i{n)\n^\v^l{n)\, 

m 



(48) 



where we remind the reader that ^fi{r) is given by 
Eq. 122). 

The result in Eq. (|48|l may be also obtained by com- 
puting the energy E{ri) of a vortex (corresponding to the 
kinetic energy of the superfluid) at position r^ in an in- 
homogeneous superfluid. Balancing the associated force 
against the Magnus "force"— gives the equation of mo- 
tion of the vortex: 



2'kTiPsZ, X Ti — V-E(r,;) = 0, 



(49) 



where the curl-free part of the vortex flow V0q is now 
regarded as part of the vortex and therefore the external 
superflow v^ in Eq. H49|l is taken to be zero. Starting 
with the kinetic energy, Eq. I|12|) (at a; = 0), expressing 
it in terms of 9y and 9a , and using the equation of motion 
for 9a, Eq. ^^, wc find 



E{r^) 



2m 



d\[ps{r){V9yf ~ p,{r){V9a) 



(50) 



with 9a the solution to Eq. (|25|l and V6'i, given by 
Eq. (|24|) . Because the curl- free superfiuid velocity cor- 
rection — 'V9a arises from the condensate inhomogeneity, 
its contribution to the vortex kinetic energy is subdom- 
inant, vanishing as 1/R^. In the dominant (first) term 
the superfluid velocity —'V9y diverges at r^, allowing us 
to ignore the r dependence of the smoothly varying den- 
sity Psir) and approximate it by its value at the vortex, 
Ps{r) ~ Ps{ri). This leads to 



t2 p 

Eiri) ~ -—psin)2n\n—, 
2m, ^ 



(51) 



which is correct to leading order in the small parameter 
^/R. When inserted into Eq. H49|) . we have 



z X Vp{ri)ln—, 

m ^ 



(52) 



in agreement (up to weak logarithmic corrections) with 
Eq. igHI) found above. 

For an axially symmetric trap, with V/i = rdrp{r), 
the equation of motion Eq. H52|l can be easily solved and 
(as advertised) gives vortex precession about the center 
of the trap: 



rt{t) = rj(xcoscjpt + ysinwpi), 



(53) 



with the direction of rotation in the same sense as vor- 
tex circulation (i.e., counterclockwise here) and with the 
angular frequency ujp given by: 



drp{n)\n—, 

mri K 

mn 2ps{ri) R 



(54a) 
(54b) 



For simplicity, in Eq. (|53|l we have chosen a particu- 
lar initial condition ri(0) = xrj. For the experimen- 
tally relevant case of a TF condensate proflle Ps{r) = 
po{l — r^/R'^), the precession frequency Eq. (|54b|l reduces 
to 



1 



i?2_ 



1 ^ 



(55) 



which agrees qualitatively with JILA experiments that 
observe such precession at angular frequency ~ h/mR^ 
that is roughly independent of 7'^. 



B. Precession of multiple vortices 

The above analysis can be easily extended to the dy- 
namics of many vortices, where each vortex moves with 
the local superfluid velocity that is due to the sum of its 
own curl- free superflow — V^^ (generated by the conden- 
sate inhomogeneity) and the superflow of all other vor- 
tices. For the simplest case of two vortices the equations 
of motion are given by 

i-i = -V0„,i(ri) + -V02(ri;r2), (56a) 

m m 

1-2 = -V0„,2(r2) + -V0i(r2;ri), (56b) 
m m 

where V0a_i(ri) is the analytic phase gradient induced by 
vortex i at position r.^ and 'V9i{rj;r.i) is the total phase 
gradient induced at position Vj due to vortex i which is lo- 
cated at position r^ . Equations H56a|) and (j56b|) are quite 
general. We now restrict attention to the case of a TF 
condensate profile and focus on the limit r^ ^ _R so that 
the rotation frequency Eq. (|55|) is approximately indepen- 
dent of radius. As we showed in Sec. lIIII 'V9a is almost 
always much smaller than V^t, (especially when r^ <ti R) 
so that the second terms in Eq. H56a|) and Eq. Ij56b|) can 
be well-approximated by neglecting the associated ana- 
lytic fields. Under these conditions we have 



ri 



r2 



h / 1 ^ , R z X (ri — r2) \ 

_(_zxriln-+ ,„^^ J^ ), (57a) 



m\R'^ 
h / 1 



e ■ (ri-r2)2 
R z X (r2 -ri) 



■(^2^ X 1-2 In , 
m\R'^ t, (r2 — ri) 



(57b) 



The two equations decouple when we change variables 
to relative (p = ri — r2) and center-of-charge (Re = 
(ri + r2)/2) coordinates: 



h /Ini?/^ 
mV i?2 



Z X p 



zx p 



Re 



—-^ In -z X Re 



(58a) 
(58b) 



so that the vortex pair center-of-charge rotates around 
the origin at the lower-critical frequency Qci 



Uc 



h R 

mR^ ^ ' 



(59) 
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and the two vortices orbit each other at frequency 



ft f I ^ R 2 



(60) 



which increases rapidly with decreasing vortex separation 
p, with the divergence at p ^ (as usual) cut off by the 
coherence length ^. 



V. VORTEX ARRAY IN A SPATIALLY 
INHOMOGENEOUS SUPERFLUID 

We now turn to the many-vortex problem, with the 
goal of computing the vortex spatial distribution in an 
inhomogeneous rotating superfluid. The superfluid ve- 
locity Vs = —'V6, measured in the laboratory frame, 
due to an array of N vortices is a solution of Eqs H13a() 
and Ijl3b|l from SeclH] By linearity of these equations it 
is given as the sum of the contributions from each vortex: 



diverges as l/|r — rj| near each vortex at r^. As illus- 
trated in Fig. 13 Vs (r) thus strongly deviates from rigid- 
body flow. In this regime, where a superfluid exhibits its 
locally irrotational quantum nature, the summation in 
Eq. (|61|l can no longer be replaced by an integration, and 
the minimization of E must be done directly over vortex 
positions, r^, rather than over a field Vs(r). In fact, as 
we saw in Sec. IIIBI this vortex discreteness manifests it- 
self even in a uniform but finite-size condensate through 
the lower-critical rotational velocity fici ~ {h/niR^) In -p 
below which no rotation is supported by the condensate. 
For a uniform inflnite condensate the problem was 
solved long ago by Tkachenkoii, who rigorously demon- 
strated that the solution is a hexagonal lattice character- 
ized by the vortex density n„o- Carrying out Tkachenko's 
exact analysis for a spatially- varying ps (r) is a formidable 
task. We shall instead develop an approximate contin- 
uum theory that nevertheless incorporates the essential 
vortex discreteness and which is valid for a smooth con- 
densate proflle Ps (r) , with accuracy controlled by V In ps . 



ve{r) 



Y^ z X (r-ri) 



r^f 



(61) A. Vortex lattice in a generic inhomogeneous 

superfluid density proflle 



with vortex positions r^ static in the frame of the nor- 
mal component (rotating with frequency fi relative to 
the lab frame). In the above, we have justifiably ne- 
glected the curl-free 'V9a contribution to the superflow 
of a vortex, studied in Sec. IIIII that, as can be seen from 
Eq. l|S5)l. scales as Vps/ps ^ 1/R and therefore gives 
a contribution that is subdominant (in our expansion in 
'V Ps/ Ps) in the large condensate limit. Moreover, since 
V^a is curl-free it cannot contribute to the rotation of the 
superfluid and thus is expected to be irrelevant for the 
rapid-rotation limit. The corresponding vortex density is 
given by Eq. H13b|l 



N 



n„(r) = (27r)-i V x V6i = ^ 5^{t - r,). 



(62) 



As discussed in Sec. Ill Al for high rotation rates Q. 
the vortex state is dense, and to a high accuracy we 
can neglect the discrete nature of vortices [embodied by 
Eqs. llM|) . l|^ ] and approximate the superfluid velocity 
Vg (r) and vortex density n„ (r) by arbitrary smooth func- 
tions that minimize the total energy, Eq. H12|l . As we 
showed in Sec. Ill Al within this continuum approxima- 
tion, the superfluid velocity is simply given by the rigid- 
body result Vs — riz x r and riyo — mfl/T:h, thereby 
providing a simple explanation for the high vortex lattice 
uniformity observed in experimcnts^'^iSiLSiS. Despite this 
agreement, it is of interest to understand the degree of 
accuracy and limitations of this uniform vortex distribu- 
tion (rigid-body) prediction. 

Away from this classical rapid-rotation limit^'^, vortex 
discreteness begins to matter and the rigid-body uniform 
vortex distribution solution clearly breaks down, as Vs (r) 



We have shown that to compute the vortex distribu- 
tion in an inhomogeneous condensate, it is essential to 
faithfully incorporate vortex discreteness in treating the 
sum in Eq. H61|) . For r near a vortex located at r^, the 
superflow is clearly dominated (see Fig. |31) by a diverg- 
ing contribution from the jth vortex, with other vortices 
making a subleading and smoothly varying correction to 
Vf?(r). To formalize this, we write r = r^ -I- Sr and split 
the sum in Eq. H61() into a contribution from the jth vor- 
tex plus a contribution due to all remaining vortices: 



voir J + Sr) = 



z X Sr 
Sr^ 



E 



Sr - r., 



(tj +Sr - r,j) 



(63) 



In the smooth correction to the superflow due to all other 
vortices the sum over vortex positions can be safely ap- 
proximated by an integral over a continuously distributed 
vortex density: 



— V9(r, + Sr) 
m 

v.(r) 



h zx Sr _ . ^ . ..^ 

T^+^s{rj+Sr), (64a) 

m or'' 

- fd^r'nAr')'-^^ 
m J (r — r )^ 



(64b) 



where hy (r) and Vg (r) are the vortex density and super- 
fluid velocity at position r, coarse-grained on the scale of 
the lattice spacing. Equation l|64b|l may be simplified by 
expressing the vortex coordinate r^ = r^ -|-u(r") with the 
r" forming a uniform hexagonal lattice at average density 
UyQ and u(r) being a static vortex displacement field. In 
terms of u(r). 



ny{r) ~ni,o(l- V •u(r)), 



(65) 
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so that the integral in Eq. Ij64b|) gives 

Vs(r) = fi[z X r — 2z X UL(r)]. (66) 

Thus, the coarse-grained part of the superflow is only 
sensitive to the longitudinal vortex displacement, related 
to u through 



UL(r) = 



VV 



(67a) 



/d2r'VG(r-r')V-u(r'), (67b) 



with G{r—r') the Green function for the Laplacian. Since 
(as we will see below) for the case of main interest of an 
axially-symmetric trap the optimum vortex lattice distor- 
tion is purely longitudinal, we may take u^ = u. Equa- 
tion H64a|l is a remarkable result that illustrates how vor- 
tices, each with a singular (f)/r superfluid velocity, add 
up to approximate rigid-body flow^^ with the (second) 
M-term characterizing deviations from it. 

To compute the vortex distribution n„(r), we express 
the energy Eq. (|12(l in terms of the vortex lattice dis- 
placement field u(r) and minimize it over u(r). To this 
end, we express the energy E of an array of vortices as a 
sum over contributions due to individual unit cells, with 
each cell associated with a particular vortes^S*^: 






(68) 



where the subscript i on the integral indicates that it is to 
be performed over a unit cell centered at r^ . In obtaining 
Eq. H68I) . we have completed the square in Eq. 112|l and 
discarded a constant term. Using Eq. H64a() for the phase 
gradient within the zth cell, shifting the integration in 
each cell via r ^ r -I- Tj and using the fact that p^ (r) does 
not vary appreciably over a cell (i.e. ps(r -I- r^) ~ Ps{i^i) 
inside a particular cell), we find 



2ujz X UL{r + rj)]' 



(69) 

The dominant contribution to the energy per cell comes 
from the diverging superfluid velocity at the center of a 
cell; thus the way in which we treat the cells' boundary 
is unimportant. Taking each cell to be a circle of radius 
a, set by the average vortex density ny(r) — l/vra^ and 
approximating the smoothly- varying field UL(r) by its 
value at the vortex position UL(r + r^) w Ui(ri), the 
integrals over the cell areas can be easily done, giving 



i?c.|l^p,(r,) 



Trln ■ 



1 



4cj2u^(r,)2 



TT^ Uv 



nvivi) 



(70) 



with the short-scale logarithmic divergence of the vortex 
kinetic energy as usual cut off by the vortex core size f . 
The vortex discreteness effects that we have empha- 
sized, arising from the singular nature of the phase gra- 
dient near the core of a vortex, are contained within the 



first term of the energy functional Eq. (|70|l and clearly 
vanish at high vortex density as Tr^^n^ ^ iJA The re- 
maining sum over vortex cells can be faithfully approxi- 
mated by an integral 



E- 



cfrny{r) .. ., 



(71a) 
(71b) 



which, after using Eq. H65|l . finally gives the energy func- 
tional of a vortex solid in an inhomogcneous rotating con- 
densate: 



E[uir)] 



2m 



(frpsir) 4uj^UL{r)^ 



+uj{l- V •u(r))ln 



^^Lu{l-V-u) 



(72) 



Armed with -E[u], the vortex distribution can be eas- 
ily computed by minimizing £'[u] over u(r), namely by 
solving ^ — 0. Since Eq. (|72|l depends only on the lon- 
gitudinal component^'* u^, (recall V • u(r) = V • UL(r)), 
we can equivalently vary with respect to u^. Thus, we 
obtain the nontrivial ground-state distortion of a vortex 
lattice from the naive uniform (rigid-body) state in an in- 
homogeneous condensate characterized by the uniaxially 
symmetric superfluid density /5s(r) = Ps(t): 



u(r) 



8ujps{r) 



Psir)ln 



euj{l-V -u) 



(73) 



where c = 1/e. For the case of Ps{r) largest at the center 
of the trap, the vortex distortion u(r) is clearly in the 
outward radial direction and purely longitudinal for the 
case of an axially symmetric Ps{t)- 

The nontrivial distortion u(r) arises as a result of the 
competition between the vortex kinetic energy (second 
term in Eq. 172(1 '). associated with vortex discreteness 
and the Magnus "energy" (first term in Eq. (|72|l ). For a 
nonuniform Ps(r) the former is lowered by shifting vor- 
tices out to the condensate edge (along — Vps), where 
Psir) and the associated kinetic energy cost is smaller. 
This is balanced by the Magnus force (which seeks to 
minimize the distortion u(r)) that is proportional to the 
difference between the vortex velocity fiz x r and the 
local superfluid velocity Vs(r) = ri[z x r — 2z x Ui(r)]. 
Since for a weak distortion the former energy is linear 
and the latter is quadratic in u(r), a nontrivial vortex 
lattice distortion, given in Eq. (|73|l . is always induced. 

Using our main result for u(r), Eq. (|73l) . the superfluid 
velocity can be easily computed using Eq. (|66|l . giving 



v<;(r) ~ riz X 



1 c 

— (ln-2-)Vlnp, 



(74) 



where for simplicity we made an approximation V • u w 
inside the logarithm of Eq. H73|l . Since typically Ps{r) 
is concave, largest at the center of trap and decreasing 
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monotonically with radius (although by tailoring a trap 
potential it can be made interestingly nonmonotonic; see 
Sec. IV B 3l belowl ■ the superfluid velocity Eq. (|7^ is in 
fact lower than the rigid-body result by an amount that 
generically increases with radius. Thus, in the rotating 
(vortex lattice) frame the superfluid velocity is in the 
direction opposite to that of the imposed rotation (see 
Fig. 13 . 

Furthermore we note that, as illustrated in Fig. |21 
Vs(r) in Eq. (|74|l exhibits radial shearing, namely the 
superfluid rotates with an r-dependent angular velocity 
and not simply as a rigid body with Vs (r) = fi'z x r and 
n' < ri (except in the case of Gaussian ps (r) , see below) . 

A conventional fluid exhibiting such a radial shear past 
a conventional (e.g., colloidal) crystal would necessarily 
exert a viscous shear stress on the crystal, thereby in- 
ducing a helical-like shear strain S^w^ (w^ an azimuthal 
lattice distortion) in the crystal. By symmetry argu- 
ments one would expect a similar helical (azimuthal 
radial-shear) distortion of the vortex lattice due to the 
shearing superfluid, in the direction opposite to the fluid 
flow. However, within our over-simplified ideal T — 
superfluid analysis (no normal, quasi-particle fluid), the 
forces on the vortex are purely radial (perpendicular to 
the superflow), and therefore only induce a radial (lon- 
gitudinal) distortion, Eq. (|751l . Whether this intrigu- 
ing, symmetry-suggested helical distortion of the vortex 
lattice will emerge from a more realistic (e.g., two- fluid 
model that includes thermally excited quasiparticles) cal- 
culation remains an interesting open question. 

Using Eq. (jT^J together with Eq. (|^ (i.e., taking the 
divergence of u(r)), or equivalently using Eq. (|66|l (i.e., 
taking the curl of Vs(r)) we can compute the correspond- 
ing coarse-grained vortex density by solving the differen- 
tial equation that ny{r) satisfies: 

n.(r) = - + ^v(^V[p.(r)ln ,,-%, ])■ (75) 

The physical picture embodied by Eq. (|75() is straight- 
forward to understand. The first term is the usual 
rigid body density discussed in the Introduction, corre- 
sponding to the vortex density in the limit of an infinite 
and homogeneous superfluid. Indeed, for uniform Psir), 
Eq. H75|l is exactly solved by the rigid-body vortex den- 
sity (i.e. Uvo = w/tt), in agreement with Tkachenko's 
resultsii. The correction to rivo (second term above) 
also vanishes in the dense vortex (or high-rotation rate, 
n^^^ « 1) limilii^, as expected from the discussion in 
Sec. Ill AI Hence, consistent with experiments^i^iSiL2iS, an 
approximately uniform rigid-body rotation correspond- 
ing to a constant vortex density rivo is predicted even in 
the case of a spatially- varying superfluid density, Psir). 
Consistent with earlier observations, since the condensate 



P's i^) 



<0, 



density variation in a trap is expected to satisfy ' — j-^ 
Eq. (f75|l predicts that the coarse-grained vortex density 
is expected to be lower than the rigid body result n^o by 
an amount that vanishes for a uniform superfluid. 



For smooth variations in Psir), we can replace ny{r) 
in the logarithm by its approximate value w/tt, yielding 
our main result 

nv{r) « - + ^v(^V[p.(r)ln-^]),(76a) 



TT Stt \ps{r) 
c\/^\nps{r), 



eoj- 



(76b) 



that relates vortex density to superfiuid density, with 

1 . c 






(77) 



Because, as discussed in Sec.|Hl in the TF limit ps{r) is 
simply determined by the trap profile with 



p,{r)^i^^-V{r) + lmn\^)/g 



(78) 



Eq. H76a() gives an explicit prediction for the vortex den- 
sity distribution in a rotating, inhomogeneous superfluid. 
It is remarkable that, in the complementary lowest 
Landau level (LLL) limit, an identical relation (with the 
exception of the coefficient c, which in the LLL regime 
is a pure number, independent of ^) emerges and was 
used by Hc44 to argue (unfortunately incorrectly^, based 
on the observed approximate uniformity of the vortex 
lattice) for the Gaussian form of the condensate profile, 
Psir) in the LLL limit. In the recent work by Watan- 
abe et al^, a TF profile was found to be the optimal 
one in the LLL limit (in agreement with experiments^) 
leading these authors to assert (consistent with our ear- 
lier prediction^-) that the vortex lattice is nonuniform in 
this regime. These latter findings were also supported 
by recent numerical solutions of the problem in the LLL 



regime 



B. 



.19.46 



Application to specific superfluid density 
profiles 



We now apply our result for nv{r), Eq. l(75|) . to a va- 
riety of condensate profiles Ps{r), realizable by tailoring 
the shape of the trap potential. 



1. Gaussian psir) 

We first note that for the case of a Gaussian conden- 
sate profile Eq. (|26|l . the differential equation for nt,(r), 
Eq. H75|l . is solved by a uniform density fiy^Q given by 



nv,G = 



In- 



TT AttB? TT^^nv,G 
LO 1 C 

m ■ 



TT AttR^ i'^UJ ■ 



(79a) 
(79b) 



This solution represents a perfect hexagonal lattice, 
static in the rotating frame, with a lattice constant 
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slightly larger than the rigid-body result. The corre- 
sponding superfluid velocity is given by the rigid-body 
form 



,(r)«f]' 



with 



n' = n 1- 



z X r, 



1 , ^ 



(80) 



(81) 



slightly smaller than the imposed rotation frequency Q. 
This result, (derived within London approximation) is 
consistent with the prediction by Ho^ (discussed above) 
of a Gaussian condensate profile for a uniform vortex lat- 
tice, derived in the complementary LLL regime. Because 
of the close similarity of this vortex configuration and the 
associated superflow to those of the rigid-body result, we 
expect difficulties for a direct experimental verification of 
our prediction in Eq. H79a() for the Gaussian condensate 
profile. 



2. Inverted parabola ps (r) (harmonic trap) 

For the most experimentally relevant case of a har- 
monic trap, (which, within the TF approximation, has 



the condensate profile ps (r) ex (i? 



[c.f. Eq. (Cni)]), 



the vortex density profile Eq. (|76a|l reduces to the result 
advertised in Sec.Q?°: 



_ , , UJ 1 W 
nv[r) ~ ^5-7^2 



.2)2 In ^2^ 



TT-.- (82) 



In Fig. n we plot Eq. H82|) for experimentally realistic 
parameters H. = .860t, R = 49/iTO (top solid curve), 
ft — .5717t, R = 31/im (middle solid curve), ft = AOilt, 
R = 25/im (bottom solid curve) along with the rigid- 
body result (dashed curve) for the case of *^Rb. Here, £, 
is taken to be the TF valueiSiS 



e 



mVLtR 



(83) 



and Vtt = 52s~^. As illustrated in Fig. 1131 our prediction 
for ny{r) compares well with the recent JILA dat£^. To 
obtain this fit, we used values for f2 and R reported in 
Ref.,^ adjusting them (from fi/f^t = .89 and R — SO/im) 
within the reported error bars (±3 and il/iTTi, respec- 
tively) to find the best fit. 

Equivalently, the above result predicts a hexagonal lat- 
tice constant {a{r) — (2/\/3n.i,(r))-'^/2) that increases with 
r, as illustrated in Figs. [T^ and ll4| where in the former 
we again compare our theory with the JILA experiments 
(with the same data as in Fig. ^. 



3. Nonmonotonic condensate profiles 

It is clear that if the only variation of the conden- 
sate density is on the scale of the trap, Eq. ((7^ pre- 
dicts a generically small (though observable-) distortion 
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FIG. 13: (Color online) Plot of a/R (lattice spacing normal- 
ized to TF radius) as a function of radius using parameters 
n/ilt = .86 and 7? = 49pm. Points are data adapted from 
Ref. [3. The dashed line is the rigid-body value of a/R. 
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FIG. 14: Plot of a/R (lattice spacing normalized to TF ra- 
dius) as a function of radius using parameters Q,/Q.t = .57 and 
R = 31/i7n. The dashed line is the rigid-body value of a/R. 



of a uniform hexagonal vortex lattice that vanishes as 
V^Xnpsir) ^ l/i?2 (for r < R)^ __To_ the extent that 



observed lattices (see, e.g., Refs. 4 5,6,7 8 9") are remark- 
ably uniform, this is a required success of our theory. 
However, to more stringently test our predictions, in the 
present section we calculate vortex lattice distortions for 
nonmonotonic (but still uniaxial) condensate density pro- 
files Ps (r) that also vary on a scale smaller than the con- 
densate's overall size R. We expect that such Ps{r) can be 
experimentally realized by confining a condensate inside 
a nonmonotonic trap potential V{r) tailored by adjusting 
a combination of magnetic and optical fields. 

For simplicity, ignoring the weak variation of the vor- 
tex density inside the argument of the logarithm, i.e., 
replacing 1 — V • u by unity, the distortion away from a 
hexagonal rigid-body array (which is a hexagonal lattice 
with lattice parameter (2/A/3r7,„o)^^^) is given by 



u(r) 






(84) 



We first consider perhaps the simplest nonmonotonic 
condensate profile with a dip in its center of spatial extent 
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FIG. 15: Main; The solid line depicts the vortex density 
nv(r)/nvo as a function of radius for the condensate profile 
Eq. II85II . The dashed line plots ps{r)/po. 



0.6 




inside Eq. (jTSJ and Eq. ^^ we find 



u(r) 



8w^2 



In- 



e^-^exp(^)-l' 



(86a) 



n„(r) = h 



1 



■In- 



P0gr2/2f2 _ ^ 
Pi 



(86b) 



2 PSlQr'-/2t^ 

I Pi 

p CP0grV2^2 _ -^p 



We plot the corresponding vortex density Eq. Ij86b|l in 
Fig. El for a large slowly-rotating condensate of ®^Rb 
atoms with parameters {R = dOOfim, O — 0.2s^^, 
Pi/ Po = 0.67, and £ = 0.3R) that are slightly differ- 
ent than typical values of present-day experiments to en- 
hance visualization of the distortion. Consistent with our 
earlier qualitative discussion, fiy (r) is largest where ps (r) 
is smallest. Near r = .5R, n«(r) actually drops below its 
rigid-body value. The corresponding distortion u(r) of 
the vortex lattice is illustrated in Fig. 1161 showing vor- 
tices displacing against the local gradient in Psir). 
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FIG. 17: Plot of condensate profile ps{r)/po (dashed line) and 
vortex density n„(r)/n„o (solid line) for the oscillatory profile 
Eq. 18711 with parameters R = 49pm and £ = 6pm. 



FIG. 16: (Color online) Distorted vortex lattice (points) aris- 
ing from the condensate profile Eq. 1851 containing a depletion 
near the origin with lines to guide the eye. 



£, that we model by 



Ps{r) ^ po- pi exp(-r /2£ ) 



(85) 



with pi < po, where the finite asymptotic density po is 
strictly speaking not a constant but varies on a much 
larger scale R and is determined by an overall large-scale 
trap with R ^ £. A qualitatively similar condensate 
profile can be generated by a "Mexican hat" -like trap po- 
tential, as experimentally demonstrated by the research 
groups of Dalibard® and Ketterle^. Using such Ps{r) 



Since u(r) is constrained to vanish at the trap center 
where Vp^ = (see Eq. (I84|l 'l. to maximize the distor- 
tions associated with a spatially- varying Ps{r) it makes 
sense to consider a condensate profile that exhibits its 
most rapid variation away from this symmetry point, i.e., 
at nonzero radii. Motivated by this we consider a trap 
potential and therefore a condensate density that oscil- 
lates with r. For concreteness we model such Ps(r) with 
a Bessel function 



Psir) = Pq + PiJoir/£), 



(87) 



with £ characterizing the length-scale of the oscillations 
and po (as in the previous example) weakly r dependent 
on a much larger overall trap scale R, with R ^ £. To 
keep the overall condensate density positive we choose 
pi < 2.48po. For this condensate profile (displayed in 
Fig. 1171 dashed curve) we compute the resulting vortex 
density hy{r) and illustrate it in Fig. El (solid curve) 
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FIG. 18: (Color online) Vortex positions (measured relative 
to R) for the profile Eq. 1871 with lines to guide the eye. The 
parameters R = 49 nm and £ — 6fim. For the upper right 
quadrant, we have shown a density plot of Ps{r) with the 
light areas denoting regions of high ps. 



FIG. 20: (Color online) Vortex positions (measured relative 
to R) for the profile Eq. I|87^ with lines to guide the eye. The 
parameters R — AQpm and £ = Afim. For the upper right 
quadrant, we have shown a density plot of pa{r) with the 
light areas denoting regions of high ps . 
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FIG. 19: Plot of condensate profile psir)/ po (dashed line) and 
vortex density n„(r)/n„o (solid line) for the oscillatory profile 
Eq. 18 711 with parameters 7? = 49/im and £ = Afim. 



for parameters (TF radius R — AQ^m, trap frequency 
fit = 52s~^, £ = 6^m, pi/po = 1.6, Q/flt = .86 and with 
£, given by Eq. IJESJ) that are consistent with recent Rb 
experiments at JILA^'^. At the location of the dip in the 
condensate density (r « 0.5i?), ny{r) exhibits a corre- 
sponding maximum. The effect on the locations of the 
vortices is quite striking, as shown in Fig. 1181 the vortex 
rings near the dip are compressed together. In the upper 
right quadrant we have included a density plot of Ps{r) 
to show how regions of low ps (r) (which are darker) in- 
duce higher local vortex density. Thus, this profile nicely 
illustrates the principle that vortices migrate to regions 



where ps is small to lower the kinetic energy. 

The case of a smaller value oi £ {i = Afj,m) , such that 
Psir) has two minima before r = i? is reached, is il- 
lustrated in Figs. ^1 and [201 again demonstrating how 
(as expected) vortices congregate near radial minima of 
Ps{t), with a local density and superfluid velocity that 
actually exceed their rigid-body values. 



C. Effect of inhomogeneous vortex density on 
condensate profile 

Thus far, to calculate the average vortex density n^(r) 
we have used the TF expression for the condensate den- 
sity Ps(r), Eq- ®! that was derived in Sec.mwithin an 
approximation of a rigid-body superflow profile. How- 
ever, as we have shown in the previous section, the 
superfiow profile is itself modified inside an inhomoge- 
neous condensate. Hence in principle we need to de- 
termine hy{r) and /Os(r) by self-consistently solving the 
coupled equations Eq. ^ and Eq. (|751) . together with 
^vi^) ^ 2^^ ^ ^s- However, as we show below, for 
a smooth trap potential, to a good approximation, the 
effect of the inhomogeneity of the vortex distribution on 
the equilibrium Ps(r) is negligible, with corrections van- 
ishing as a higher power in Vps/ps ~ 1/-R- 

To show this, we repeat the steps of Sec. lV Al but now 
using the full bosonic energy density Eq. O instead of 
keeping only the 6'-dependent terms (i.e. as in Eq. H68|) ). 
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This gives: 



E^-ld^rpAr) 



In 



TT^^Hy 



^Lo^ULirY' 



(fr 



[V - \mn\^ - /i)p.(r) + \ps{vf] ,(88) 



where we remind the reader that n^ is related to u^ via 
Eq. (jHHJ, and, in the spirit of the TF approximation, 
we have still neglected the subdominant (away from the 
condensate edge, r <C -R) Vp^ contributions to the en- 
ergy. The corresponding, coupled Euler-Lagrange equa- 
tions {SE/Sul — 0, 6E/6ps = 0) that determine n„(r) 
(through ul) and /9s(i") are given by: 



now study long-wavelength vortex fluctuations about this 
slightly inhomogeneous equilibrium vortex state. To this 
end we derive a vortex lattice elastic energy by expanding 
the total energy Eq. (|72|) to harmonic order in small dis- 
tortions e(r) about the ground-state configuration Uo(r) 
[satisfying Eq. (|75|) ]. defined by u(r) = Uo(r) +e(r). We 
find 



Eei[eir)]=^l d'rp^ir) 



2m 



(frps{r) 



^^'ei-^ 



uj (V-e) 



Au\^ 



2 1-V-uo 



,(91a) 



^__(V.e)- , (91b) 



UL(r) 



1 '^/'-W.jj^. c 



%uj ps{r) ^'^uj' 



(89a) 



Ps{r) ^-[p-bhn- (Vir) - lmn\^)] 
9 2 ' 

- i \2mn^ul - nnV ■ u in ^1 , (89b) 
9 C^w 



where, as before, in Eq. H89a|) we have replaced 1 — V-u ~ 
1 in the argument of the logarithm of Eq. (|73|l . and de- 
fined a parameter 6 = i In 7^ . While the equation for u^ 
remains unchanged, the equation for Psi'r) is modified by 
the vortex kinetic energy, that, as expected, suppresses 
it. 

Although these coupled equations are in general quite 
nontrivial, for a smooth trap the distortion u(r) is small, 
vanishing with ^ps/ps ~ 1/R. Hence, they can be 
straightforwardly solved by iteration in this small param- 
eter. To lowest order we can simply neglect the vortex 
lattice distortion u(r) inside the /Os(r) equation, obtain- 
ing our previous symmetric TF result 



where in Eq. (|91b|l we have approximated uq ~ 0. 

This compressional energy must be augmented by 
the elastic energy of the vortex lattice shear deforma- 
tion, that is clearly missed at this level of a coarse- 
grained (density-functional) approximation. We follow 
Baym and ChandleeiS and fix the shear modulus using 
Tkachenko'siii^ exact result for the uniform condensate, 
that we expect (given the high uniformity of the vortex 
lattice observed in experiments and demonstrated here 
analytically) to be a good approximation even for our 
case of a spatially- varying Ps(j). This yields 



Ee 



2m 



cPrps{r) 



Aujhl 



-{V-ef 






.(92) 



Ps{r) 



9 



bhn - -m{n^t 



n^y 



(90) 



corrected only by bhn^ that is small compared to fi. 
Higher order corrections to Psir) are easily obtained by 
an iterative proceedure, thereby generating a perturba- 
tive expansion in 1/R. Hence, in the high rotation limit, 
n ^ fit (corresponding to a divergent i?(ri), Eq. (|11|) ') 
the TF condensate density profile Eq. @ is an accurate 
description, in agreement with recent results obtained 
within the lowest Landau level approximationiSiiS. 



The dynamics of long-wavelength elastic vortex fluctu- 
ations is governed by the balance of the Magnus "force" 
against the "elastic" force associated with the distortion 
e(r, t) of the vortex lattice 



SE,, 



Lu 6e(r,t) 



27rhpsZ X e(r, f) = 0, 



(93) 



VI. ELASTIC THEORY AND TKACHENKO 
MODES 

Having determined the equilibrium vortex distribution 
in a spatially inhomogeneous rotating condensate, we 



where 



_ de 

dt 



and the inverse of the average vortex den- 



sity factor l/uyo — tt/w arises from the Jacobian relating 
the variation with respect to r^ to the functional varia- 
tion with respect to ^(ji). Carrying out the functional 
derivative yields the following equation for e(r,i): 
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PsZ X 



m J am 



r 



A full solution of Eq. (|94|l is beyond the scope of this 
paper and we leave it for future research. In spite of its 
complexity it is reassuring that for a uniform condensate 
(r-independent ps) the divergence and curl of Eq. (|94|) are 
equivalent to Eqs. (71) and (72) of Ref.|42, respectively. 
In the following, we proceed by making a simple analytic 
approximation to Eq. (|94|) to bring out the leading order 
effect of the condensate inhomogeneity on the Tkachenko 
waves. 

For long- wavelength Tkachenko modes, with - 3> fc ^ 
-^, it is clear that the first term on the right side of 
Eq. (|94ll (which is purely longitudinal) dominates over 
the second since ui ex 1/a^. Taking the divergence of 
both sides of Eq. (|94|l and neglecting the subdominant 
second term, we find 



V xe= -2nV 



(95) 



to leading order in a double expansion in spatial deriva- 
tives of e and ps- The transverse part of Eq. H94|) is 
obtained by taking its curl, which leads to 

V X [psi X e] = -^ [V X (psV^e) 
Hm L 

+V X ([Vps ■ V]e + [z • {Vps X V)]z x e) . (96) 

Henceforth, for simplicity, we restrict our attention to 
the case of a Gaussian Ps{r), Eq. (|26|l . which satisfies 
Vps{r) = —vps{r)/E?. To leading order in spatial gra- 
dients, the quantity in parentheses in the second line of 
Eq. l|^ reduces to 



[Vps ■ V]e + [z • (Vp, X V)]z X e 






(97) 



Inserting Eq. (|97|) into Eq. H9()|l and expanding the result 
to leading order in derivatives of ps , we have 

V.ez^--A^Vxe-Av2(Vxe), (98) 

where we have used V x W^e = V'^(V x e). Combining 
this result with Eq. (|95|l . we finally obtain an equation 
for the transverse vortex fluctuations V x e: 



dt^ ^ '~ ^m\ i?2 



V" (V X e). (99) 



The spatial and temporal Fourier transform of this lin- 
ear wave equation then immediately gives the Tkachenko 
mode dispersion ^^(fc): 



Vlrik) = ^Jh^/Am^jB - R- 



(100) 



that, in the uniform limit i? — > cx), recovers the stan- 
dard linear dispersion result ^T{k) — ky/hfl/Am (see 
Refs. 4^481). Given the small gradient in Ps{r) expan- 
sion that led to this result, the wavevector is obviously 
limited to the physically sensible range 1/a > k > 1/R. 
No doubt a more detailed analysis of Eq. (|M|) is needed to 
obtain a more accurate prediction for Tkachenko eigen- 
modes and their dispersion in an inhomogeneous conden- 
sate characterized by a generic condensate profile Ps{r)- 
We leave such analysis, as well as an extension of Eq. I|94|) 
that includes condensate density fluctuations, as dis- 
cussed in the recent work by Baym^ (see also Ref. 50), 
to future research. 



VII. SUMMARY AND CONCLUSIONS 

To summarize, we have presented a London theory of a 
vortex state in an inhomogeneous rotating Bose-Einstein 
condensate. Our most important result (Eq. O) is 
the relation between the vortex density fit,(r) and the 
superfluid density Psif). When applied to a harmonic 
trap, this result provides a simple explanation for the ob- 
served highly regular vortex arrays in such strongly inho- 
mogeneous condensates, and the observed Thomas-Fcrmi 
parabolic condensate profile. This relation also predicts 
a slight inhomogeneity in the vortex spatial distribution, 
that has recently (since our prediction) been observed 
in experiments and simulations from JILA— . As we dis- 
cussed in the main text, this relation between n„(r) and 
Psir) can be more stringently tested by studying vortex 
lattice distortions induced by nonmonotonic condensate 
profiles, tailored with various trapping potentials. 

As an important digression, we also studied the su- 
perflow in the vicinity of an isolated vortex, calculating 
the superfluid velocity distortion induced by a spatially- 
varying condensate. The associated additional contri- 
bution to the superflow is directed orthogonal to the 
displacement of the vortex away from the trap center, 
thereby providing a simple explanation for precession 
of such vortices about the trap center with radially- 
independent frequency. 

With respect to our vortex lattice result, it is impor- 
tant to stress that our prediction of the deviation of 
the vortex density from the spatially independent (rigid- 
body) result is more than just a suppression of the to- 
tal vortex number (a somewhat trivial and experimen- 
tally difficult ^^ effect to detect, that follows directly from 
the existence of J7ci, as can be seen from arguments in 
Sec. Ill B1 and Fig. Q. In particular, it is not a simple 
surface effect, such as, e.g., a missing vortex ring at the 
edge of the condensate, as for example observed in simu- 
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lations of Feder and Clark^^ and predicted theoretically 
for uniform supcrfluids some time ago^. 

Our prediction of a radially increasing vortex lattice 
spacing across the bulk of a condensate (that does not re- 
quire a precise measure of fi^^), shows quantitative agree- 
ment with recent JILA data^, exhibiting the same ~ 2% 
lattice distortion. 

We conclude by noting that a number of important ex- 
tensions of our work remain. Based on symmetry consid- 
erations we have suggested that an azimuthally directed 
vortex lattice radial-shear distortion drUij, will be induced 
by the radially-shearing superfluid. Clearly, it is impor- 
tant to determine theoretically whether such interesting 
chiral vortex lattice distortion is indeed present. This 
should be possible to assess within either a two-fluid hy- 
drodynamic model or by incorporating thermally-excited 
Bogoliubov quasi-particles into our theory. A more care- 
ful treatment of vortex discreteness (that goes beyond 
our simple density- functional-like approximation), from 
which a finite vortex lattice shear modulus must emerge is 
also highly desirable. This would allow a first-principles 
treatment of Tkachenko modes, that is more precise than 
that presented here. Finally, the extension of our work 
to three dimensions, where a nontrivial z-dependence of 
the vortex-line density profile is expected, is another im- 
portant future direction. 
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by an exact calculation of vortex precession in a homoge- 
neous condensate confined to a finite "bucket" of radius 
R (see Ref. 123) . To model this system we take ps be a 
constant po for r < i? and zero for r > R. For this case, 
the single- vortex problem is solved by the method of im- 
ages to enforce the boundary condition of no outward 
current (i.e. V6' must be orthogonal to the boundary)^. 

It is straightforward to verify that a single vortex at 
Tq in the bucket geometry has the solution 



V6'(r) = 



z X (r-ro) z X (r - r^) 



(r-ro)2 



(r- 



(Al) 



where r,; = roR^/rl is the "image charge". Thus, al- 
though the condensate is homogeneous in the bulk, the 
nontrivial boundary condition necessarily introduces a 
curl-free (in the bulk) correction V^a to the superflow 
(second term in Eq. ljAl|l 'l. that has a similar physical 
effect as that due to the bulk condensate inhomogeneity 
studied in Sec. IIVI 



As for an inhomogeneous condensate, the background 
superflow at the vortex and thus the vortex dynamics are 
entirely determined by the image charge term. Evaluat- 
ing the associated contribution at Tq and using Eq. 147|) 
we find 
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APPENDIX A: VORTEX PRECESSION IN A 
UNIFORM SUPERFLUID 



Our derivation in Sec. IIVI of vortex precession in an 
inhomogeneous trapped superfluid can be complemented 



To 



R^-r^ 



(A2) 



showing that a vortex located off-center in a finite 
"bucket" will precess about the origin in a sense simi- 
lar to that embodied in Eq. l|^ and Eq. l(5^ . with the 
frequency near the center of order Ti/raR^, correspond- 
ing to a unit of angular momentum for a boson at the 
boundary of the condensate, as found for an inhomoge- 
neous condensate. 
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